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Spin-only descriptions of the half-filled one-band Hubbard model are relevant for a wide range 
of Mott insulators. In addition to the usual Heisenberg exchange, many new types of interactions, 
including ring exchange, appear in the effective Hamiltonian in the intermediate coupling regime. 
In order to improve on the quantitative description of magnetic excitations in the insulating anti- 
ferromagnetic phase of copper-oxide (cuprate) materials, and to be consistent with band structure 
calculations and photoemission experiments on these systems, we include second and third neighbor 
hopping parameters, t' and t" , into the Hubbard Hamiltonian. A unitary transformation method is 
used to find systematically the effective Hamiltonian and any operator in the spin-only representa- 
tion. The results include all closed, four hop electronic pathways in the canonical transformation. 
The method generates many ring exchange terms that play an important role in the comparison with 
experiments on La2Cu04. Performing a spin wave analysis, we calculate the magnon dispersion as 
a function of U,t,t' and t" . The four parameters are estimated by fitting the magnon dispersion to 
the experimental results of Coldea et al. [Phys. Rev. Lett. 86, 5377, 2001] for La2Cu04. The ring 
exchange terms are found essential, in particular to determine the relative sign of t' and t" , with the 
values found in good agreement with independent theoretical and experimental estimates for other 
members of the cuprate family. The zero temperature sublattice magnetization is calculated using 
these parameters and also found to be in good agreement with the experimental value estimated 
by Lee et al. [Phys. Rev. B 60, 3643 (1999)]. We find a value of the interaction strength t/ ~ 8t 
consistent with Mott insulating behavior. 



I. INTRODUCTION 

High temperature superconductors have challenged al- 
most every traditional concept and method of condensed 
matter theory. Ill understood issues concerning Fermi 
liquids and quantum critical behavior for example, may 
be intimately related to the problem of superconductivity 
itself. From the very beginning, it has therefore appeared 
important to develop a better understanding of quantum 
magnetism in general, given that the parent compounds 
of high-temperature superconductors are insulating anti- 
ferromagnets. For example, although these compounds 
are Neel ordered, it was proposed early on that, in two 
dimensions, the ground state of the Heisenberg antiferro- 
magnet might be very close to being a resonating- valence 
bond state rather than a Neel ordered statoi. This turns 
out not to be true, but it took some time to establish 
numerically that the ground state is indeed ordered^i^, 
in agreement with experimental^. 

Issues relating to the quantum magnetism of the 
copper-oxide materials (cuprates) have resurfaced in the 
last few years. Detailed neutron scattering experi- 
ments on the parent high-temperature superconductor 
La2Cu04 illustrate that the question of which model best 
describes these compounds^ can now be addressed with 
precision. For example, while the magnetic properties of 
Copper Deuteroformate Tetradeuterate can be described 



in detail using the Heisenberg model^,, La2Cu04 exhibits 
clear deviations^ from it. 

Specifically, the observed spin wave dispersion along 
the antiferromagnetic Brillouin zone boundary is not pre- 
dicted by the Heisenberg model. One can in principle 
account for these differences by introducing phenomeno- 
logical two-spin and multi-spin exchange constants into 
the spin Hamiltonian. These spin Hamiltonians can be 
deduced, as shown by Hhia.c^'^^, from general consider- 
ations of symmetry and permutation operators and the 
value of the exchange constants can be fixed by fitting to 
the experimental data. However, in our opinion, a more 
interesting approach is that taken by Coldea et alM-. The 
authors of Ref. take the point of view, which we share, 
that the one-band Hubbard model is a more fundamen- 
tal starting point for the description of the copper oxy- 
gen planes in these materials if one wants to connect the 
magnetism at the microscopic level with electronic cor- 
relation effects. Although in this paper we concentrate 
on half filling, it may also be valid on moving away from 
half-filling, into the region where the fermions eventually 
pair to give high-temperature superconductivity. 

The one band Hubbard Hamiltonian iJjj is one of 
the simplest lattice models of interacting electron that 
admits a "rich" phenomenology resembling that of the 
cuprates. A key feature of this model is that it does de- 
scribe both the insulating and metallic phases and the 



(Mott) transition between them. In interaction-induced 
insulators, the band is half-filled and insulating behav- 
ior occurs because of interactions, not because of band 
folding induced by long-range order. However, in two- 
dimensions a half-filled band can be insulating because of 
a pseudogap induced by fluctuating antiferromagnetism, 
or because of the Mott phenomenon. In the former case 
(weak coupling), the pseudogap appears when the ther- 
mal de Broglie wavelength ^th is less than the antifer- 
romagnetic correlation lengtbii ^, while in the latter 
Mott insulating case (strong coupling), the system is in- 
sulating even when this condition is not satisfied. In 
La2Cu04 neutron measurements have been done only in 
the regime where ^ > ^th, so that antiferromagnetic cor- 
relations could be the source of the insulating behavior. 
In a very recent paper— Comanac et al. used single- 
site dynamical mean-field theory (DMFT) to calculate 
the optical conductivity of La2Cu04. The authors of 
Ref. found that La2Cu04, or at least its optical con- 
ductivity, is best parametrized by a [//t smaller than that 
necessary for the Mott transition in single-site DMFT. 
Consequently, the authors of Ref. [l^l argue that corre- 
lations in the cuprates may not be as strong as gener- 
ally believed and that these materials may not be Mott 
insulators but rather, that the development of antiferro- 
magnetism is necessary to drive the the insulating state. 
On the other hand, c/wsier-DMFT calculations, which 
are more accurate in two dimensions, find that the inclu- 
sion of short-range antiferromagnetic correlations, which 
is possible in multi-site calculations, reduce considerably 
the critical (f//i)Mott necessary for the Mott transition 
compared with single-site OMFTi^iiiii^O^. 

This leads to one of the key questions in this field. Do 
the strange properties observed in underdoped cuprates 
at finite temerature arise from interaction-induced local- 
isation (Mott insulator) or from competing phases? If 
interactions are not strong enough to lead to a Mott in- 
sulator, competing phases arising at weak coupling are 
the only remaining possibility. 

To find parameters appropriate for high-temperature 
superconductors, Coldea et al. used earlier theoretical re- 
sultsiiii^ relating the parameters of the Hubbard model, 
hopping energy t and on-site interaction [/, to the ex- 
change constants in an effective spin-Hamiltonian. The 
parameters in the Hubbard model were deduced by fit- 
ting the resulting effective low-energy spin theory to the 
spin-wave dispersion relation. This is quite a remarkable 
result. Indeed, had the Heisenberg limit accurately de- 
scribed the Hubbard model, only the value oit^/U could 
have been deduced from neutron experiments. The ap- 
pearance of corrections to the Heisenberg model to higher 
order in powers of t/U allows one to obtain t and U 
separately, the former being an effective electronic band 
quantity. 

However, to be consistent with a general message pro- 
vided by Angle Resolved Photoemission Spectroscopy 
(ARPES)^i^i2ii^ experiments, optical experiments^^ 
and with results from band structure calculation a^^'^^ on 
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FIG. 1: (Color online) Hopping processes and resulting spin 
interactions, a) shows the different hopping processes char- 
acterized by parameters t, t' and t" in the Hubbard model 
considered in this paper, t, t' and t" are hopping param- 
eters between first, second and third nearest neighbors, re- 
spectively, b) At half-filling, and when t' — t" = 0, canonical 
perturbation theory leads to order t*' /U'^ to an efi^ective spin- 
1/2 Hamiltonian, , characterized by first ( Ji), second ( J2) 
and third (J3) nearest neighbor exchange interactions as well 
as a 4-spin ring (cyclic) exchange interaction with coupling 
strength Jc- c) In a large S expansion, can be recast 

to order 1/S as an effective spin Hamiltonian which only in- 
volves bilinear (pairwise) spin-spin exchange interactions be- 
tween first ( Ji*^), second {J2^) and third {Jz^) nearest neigh- 
bors. To order 1/S, the effect of the ring exchange of strength 
Jc merely renormalizes the first and second nearest-neighbor 
exchange as discussed in Ref. . d) illustrates an example of 
a four hops (ring exchange) electronic process that involves t' 
and t" and contribute to order t^t't" /U'^ to the spin Hamilto- 
nian Hs*\ e) shows that the ring exchange term illustrated 
in d) introduces in the 1/S approximation of H^"^ a fourth 
nearest neighbor exchange, as well as renormalize the 
Jf, Jf and Jf of c). f) illustrates all the additional 
exchanges beyond those shown in c) and which are generated 
by all the hopping processes involving an allowed combination 
olt, t' and t" to order 1/U^. 



a wide variety of quasi two-dimensional cupratc materi- 
als, it is necessary to include in the Hubbard model both 
second {t') and third {t") nearest- neighbor hopping, as 
they have been found to have sizable values in all the 
above studies (see Fig. [T]). In this paper we thus ad- 
dress the following questions: What is the effect of these 
additional hopping constants on the magnetic properties 
of the cuprates, and on La2Cu04 in particular? Can 
their "effective" values and that of U within the Nccl or- 
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der phase be determined by comparing experimental spin 
wave dispersion data with theoretical calculations? 

As a first step to address the question of the strength 
of t, t', t" and U in La2Cu04, we incorporate the effects 
of these extra hopping constants into a one-band Hub- 
bard model from which we derive an effective spin-only 
Hamiltonian description of La2Cu04. As we shall see, 
the effects of t' and t" are rather subtle and can even 
sometimes compensate each other. It is only when these 
hopping terms have been included in the derivation of 
the effective spin Hamiltonian, with all the subsequent 
new ring exchanges, that one can extract more meaning- 
ful values oi t — t' — t" — U from comparisons of theory 
with inelastic neutron scattering experiments. In partic- 
ular, the values of t and U obtained by Coldea et al^ 
place the underlying Hubbard model in the intermedi- 
ate coupling regime where the band width and interac- 
tion energies are comparable and hence in the regime 
where one might expect an insulator to metal transition. 
Our more detailed analysis leads to a larger values of [/, 
repositioning material relevant model parameters more 
clearly within the Mott insulating regime, as we will dis- 
cuss more thouroughly at the end of this paper. 

Although we analyze experiments on La2Cu04, pro- 
viding a specific materials context to our study, we be- 
lieve that the parameters of the one-band Hubbard model 
that we extract below should be characteristic of the two- 
dimensional copper-oxygen planes of cuprate supercon- 
ductors in general: specifically, band structure calcula- 
tions do reveal variations of band parameters from one 
compound to the nex t^^i^^ , but overall, the band struc- 
ture of Cu02 planes is quite similar from one compound 
to the next^. Even away from half filling, whether for 
hole or electron doped materials, the one-band Hubbard 
model for Cu02 planes seems to contain much of the 
Physics of high-temperature superconductor o^^'^^ . Our 
calculations should therefore have wider ranging applica- 
tions than the detailed comparison with experiments on 
pure La2Cu04 presented below. 

In the next two subsections of this Introduction, we 
describe in more detail the central issues and key results 
of this paper. There are two main parts to our work. In 
the first part. Section [III we derive the effective Hamil- 
tonian to order l/U^^ including t ~ t' — t" . This deriva- 
tion is exact and completely general. In the second part. 
Section III, we determine the spin wave energies, taking 
into account the quantum correction to the classical fre- 
quencies through the spin wave velocity renormalization 
factor, Zc{k), to lowest order in 1/S. We then determine 
a parameter set {t,t' ,t" ,11} by fitting the experimen- 
tal magnon excitation spectrum of Rcf. Q. The result- 
ing set is in good agreement with values from ARPES 
experiments and band structure calculations on various 
cuprates. Further, we calculate the zero temperature 
staggered magnetization order parameter using the set 
{t, i', t", U} and compare it with experimental values. A 
more complete solution of the spin dynamics for the effec- 
tive spin Hamiltonian derived herein, including magnon- 



magnon interactions, would most likely require the use 
of numerical techniques^, since it would be an arduous 
task to include in the self-energy corrections the multi- 
tude of magnon creation and annihilation terms present 
beyond the non-interacting approximation. 



A. Spin-only description of Hubbard model and 
La2Cu04 

In many strongly correlated quantum mechanical sys- 
tems, the separation of energy scales allows the develop- 
ment of a low energy effective theory, through the integra- 
tion over the high energy degrees of freedom. In this con- 
text, spin-only descriptions of strongly-correlated elec- 
tron systems are an excellent example. Here, in the limit 
of strong electron-electron interaction, elimination of the 
states with one or more doubly occupied sites reduces 
exponentially the dimension of the relevant low energy 
Hilbert space, facilitating the theoretical and numerical 
description of the problem. Through this procedure, the 
regime of particular interest, where energy scales actu- 
ally compete, can be approached perturbatively, via the 
introduction of a small parameter, the ratio of the kinetic 
to potential energy scales. However, the price one pays 
for this dimensional reduction is the generation of longer 
range and many particle interactions, as one moves into 
this intermediate regime. The derivation of such a the- 
ory can be achieved using different methods, leading to 
distinct expressions for the effective Hamiltonian, which 
can be shown to be equivalent through the application of 
a unitary transformation^. One of these methods, the 
canonical transformation (CT) ^^i'^^i^'^ , which we adopt 
below, is a convenient and systematic way of expanding 
a model like the Hubbard model to any order in the per- 
turbation parameter. 

In the regime where t/U remains a small parame- 
ter suitable for a perturbation theory, but where one 
moves away from the strictly Hcisenbcrg limit {t/U — > 0), 
the pairwise exchange interactions are joined by ring, or 
cyclic exchange terms that couple more than 2 spinsi^. 
Large enough ring exchange terms have been found to 
drive some model systems into exotic, intrinsically quan- 
tum mechanical ground state a'^^'^^i'^^ . Indeed, there is 
currently a rapidly growing interest in the study of ef- 
fective many-body non pairwise spin— spin interactions 
which lead to non-trivially correlated states in quantum 
spin systems— 

In most theoretical work, the effort has so far focused 
on models with only the nearest-neighbor Heisenberg ex- 
change and ring exchange. However, when starting from 
a microscopic fcrmionic Hamiltonian the strength of the 
ring exchange coupling is determined from the small pa- 
rameter t/U and is not a free parameter, nor is it the 
sole higher order term arising in the resulting spin only 
theory. To the same order in perturbation, t{t/U)^, the 
canonical transformation also generates second, J2, and 
third, J3, neighbor spin-spin interactions through pro- 
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cesses involving four electronic hops. The ring exchange, 
Jc, is thus only one of a set of spin interactions generated 
to this order— i^i^ and all should be taken into account. 
We remark that in the general approach of Dirao^i^, the 
coupling constants multiply spin permutation operators, 
rather than spin operators themselves and thus have a 
different definition from those here. One finds that per- 
mutation of four spins lead to both four spin ring ex- 
change terms and pure two-spin exchange when written 
in terms of spin operators^^. 

Within the t — U Hubbard model, the inclusion of vir- 
tual hopping pathways of n hops thus introduces fur- 
ther neighbor spin couplings of order t{t/U)'"'~^ . For a 
compound with a band structure characterized by sev- 
eral tight-binding parameters of comparable value, one 
should consider the role of direct further neighbor hops 
over length scales comparable with the nearest neighbor 
pathways of length n, as well as the multiple hop terms. 
Hence, to develop a theory accurate to order [l/UY' one 
should include direct hopping to first, second and third 
nearest neighbors, characterized by energy scales t, t' and 
i", respectively. These hopping processes generate sec- 
ond and third neighbor Heisenberg exchange terms that 
compete with those generated by the four hop processes 
discussed above. There is evidence that this is the sit- 
uation for the Cu02 planes of parent superconducting 
materials: for example, in Sr2Cu02Cl2, which has a half 
filled band, the hopping constants t' and t" are estimated 
to be t' /t ^ —0.3, t" /t ~ 0.2, from comparisons of exact 
diagonalization with photoemission experiments^. Very 
similar values are found from comparisons with other 
types of calculations^ and from band structure calcu- 
lations for YBa2Cu0307i^. In fact these values are typi- 
cal for insulating materials with copper-oxygen planes, 
including both hole and electron doped cuprates^Li^. 
For La2Cu04, band structure calculations^!^ suggest 
smaller values of t' /t and t" /t. ARPES experiments 
are not available for this compound but they are for 
the doped system La2-a;Sra;Cu04 that becomes a high 
temperature superconductor. According to these exper- 
imentsi^i2£i^i2^, t' /t and t" /t deduced from the shape 
of the Fermi surface are close to band structure values, 
although they are larger by about 30% for very lightly 
doped compoundi^. We will return to this issue later. 

Given this range of values for t' and t" one could ex- 
pect the resulting exchange terms of order t'{t' /U) or 
t"{t" /U) to be of similar magnitude to those of order 
t{t/UY . This is one of the motivations for the work pre- 
sented in this paper where we introduce all terms in a 
spin-only Hamiltonian that depend on i, t' and t" and 
that arc generated up to order . We assess their im- 
portance, making particular reference to the magnetic 
properties of La2Cu04. We show that the canonical 
transformation methocti^i^i^ can be easily generalized 
for the t — t' — t" — U Hubbard model, giving a complete 
spin-only description for arbitrary first, second and third 
neighbor hopping, correct to order IjU^. 



B. Spin- wave analysis and comparison with 
experiment. 

It has been known for some time that a spin wave 
analysis of the spin 1/2 Heisenberg antiferromagnet on a 
square lattice, taken to leading order in l/S, reproduces 
the zero point quantum fluctuations of the staggered mo- 
ment to a good approximation^. This is because the 
second order terms renormalize the classical magnon fre- 
quencies but do not induce magnon-magnon interactions 
and make zero contribution to the magnetization. In 
this paper, we therefore take it as a reasonable first ap- 
proximation to calculate quantum spin fluctuations to 
lowest order in l/S* and include the first order correc- 
tions to the classical frequencies. We have not addressed 
the question as to whether other physics, in particular 
magnon-magnon interaction is introduced to second or- 
der in 1/S. 

We compare our results with the inelastic neutron scat- 
tering data of Coldea et al^ and with the spin wave calcu- 
lation therein where only nearest neighbor hopping t was 
considered. Wc find that including further neighbor hops 
allows for a more convincing description of the neutron 
data, with the set of parameters {t,t' ,t" ,U) that best 
describe the data being comparable to other estimates 
for cuprate materials. Our analysis also gives a new es- 
timate for the total staggered magnetic moment which 
is in good quantitative agreement with the experimental 
estimate^. 

We now describe in more detail the main features of 
the data that must be explained and how going beyond 
the nearest-neighbor hopping provides a better model. 
The inelastic neutron scattering data of Coldea et al^ 
provide an extensive description of the magnon disper- 
sion in La2Cu04 over the whole of the 1'^' Brillouin zone. 
The analysis of Ref. clearly shows that the nearest 
neighbor spin 5* = 1/2 Heisenberg Hamiltonian docs not 
reproduce the magnon dispersion over the whole zone. 
Most notably, while a spin wave analysis of the Heisen- 
berg model to order 1/S' gives a flat magnon dispersion 
over the interval qbz G [(""jO), (7r/2,7r/2)] (along the an- 
tiferromagnetic zone boundary) , the experimental results 
show negative (downward) dispersion over this region (see 
Fig.[3|). Going beyond l/S for the Heisenberg model does 
not explain this negative dispersion feature^^^. Rather, 
introducing interactions between magnons gives rise to 
a small positive (downward) dispersion over this qbz in- 
terval^i^. Including terms up to t{t/U)^, which gener- 
ates second and third neighbor and ring exchange inter- 
actions, does however provide the characteristic down- 
ward dispersion, and this is what allowed Coldea et al. 
to fit the experimental data within a, 1/S spin wave anal- 
ygjg47j48^ Recent analysis of triplon excitations in the 
copper oxide based ladder material La4SrioCu2404i also 
suggest that ring exchange interactions are present and 
are of a similar amplitude^ to those found in Ref. Q for 
La2Cu04. 

It is important for the rest of the paper to expand 
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briefly on this negative dispersion and to give an inter- 
pretation for it: the spin-only Hamiltonian, Hs'^\ ob- 
tained from the Hubbard Hamiltonian, Hu, to order 
f^/U^, contains first, second and third nearest-neighbor 
exchanges, Ji, J2 and J3, respectively, as well as a four- 
spin ring exchange Jc (see Eq. ([T6l) ). Within the Holstein- 
Primakoff spin-wave approximation, all terms in Hs'^\ 
including the ring exchange, contribute to order l/S 
quadratic magnon creation/annihilation terms to the re- 
sulting quadratic magnon Hamiltonian q^^d- 
order, the spin wave approximation thus eliminates the 
original multi-spin (ring exchange) nature of the spin 

Hamiltonian, Hg'^K As a result, H^"^^^^^ could equally 
well have been derived from a different spin-only Hamil- 
tonian, i?s^, containing solely bilinear, or pairwisc, spin 
exchange terms only. would have exchange terms of 
the form J°^S{0) ■ S'(r„), which couple a reference spin 
S{0) to an 71— nearest neighbor spin at r„ 47,48^50 ^ j,^^ 
relationship between the J^^ exchange couplings of 
and the set {Ji, J2, J3, Jc} of i/i''^ is^: 

Jf - J2 - JcS' } ■ (1) 
Jf = J3 J 

where S = 1/2, Ji = (4<VC/)(1 - 6t^/U^), J2 = J3 = 
4t^/U^ and Jc = 80f*/f/^, and where the convention 
J^^ > means antiferromagnetic exchange—. We note 
from Eq.[l]that Jc leads to a "renormalization" of J2. In 
particular, for the t~U Hubbard model, = —Idf^/U^ 
{S = 1/2) is negative for all t/U. The negative value 
of Jf'^ favors ferromagnetic correlations across the di- 
agonal of a plaquette and gives a downward dispersion 
along qez G [(""i 0), (7r/2, 7r/2)] — , as required. In the 
limit t/U ^0, one recovers the nearest- neighbor Heisen- 
berg 5=1/2 Hamiltonian and the dispersion becomes 
flat along the zone boundary^. 

Having provided an explanation for the origin of the 
negative (downward) dispersion along Qbz, we can now 
discuss the effect on the spin wave dispersion of including 
second nearest- neighbor hopping t' . When passing from 
the Hubbard model, including to a spin only model, 
the leading order effect of t' is to generate an antifer- 
romagnetic exchange between second-nearest neighbors 
that modify J2 in ifi^'': 

J2 ^ j2 + i{t')yu = uyu^ + A{t')yu (2) 

which in turn modifies in the 1/5 spin- wave approx- 
imation in Eq. p]): 

Jf ^ jf + 4{t')yu = \(t'fiu-m^iu^. (3) 

So, while for = the dispersion is downward for all tjU , 
a nonzero t' competes with this trend since the second 
nearest-neighbor exchange generated to order 4(t')^/?7, 
is antiferromagnetic and frustrating. This means that a 
sizable increase in t/U is required to achieve a good fit to 



the data of Ref. [a| when t' is included in the spin-only 
Hamiltonian to this order. The same observation has 
been made from a random phase approximation calcula- 
tions^!^ and also from quantum Monte Carlo calcula- 
tions^, both on the half-filled one band Hubbard model. 
This evolution puts the best fit value for t/U into the 
intermediate coupling regime and perhaps very close, or 
maybe even beyond the critical value for a metal insula- 
tor transition. If we estimate this critical value from the 
value of U at which a finite gap in the density of states 
persists at finite temperature, then Fig. 5 of Ref. (56j 
suggests that for t' = t" = on the square lattice, the 
critical value is t/U ^ 1/6 = 0.166. In the variational 
cluster approximation, which overestimates the effect of 
interactions, the critical value of t/U is larger but stud- 
ies as a function of \t'\ suggest that frustration leads to 
a decrease of this valued. The same trend as a func- 
tion of frustration [t') occurs on the anisotropic trian- 
gular latticoSI. Hence, with nonzero t' included in the 
model and presumably at play in the real material, it is 
ultimately important that successful fits to experiments 
on parent insulating compounds lead to smaller values of 
t/U than those where t' = 0, contrary to what is found in 
the references cited above where the third nearest neigh- 
bor hopping t" is equal to zero. 

As discussed above, the third neighbor hopping term t" 
is estimated to be of the same order as t' in a number of 
cuprate materials and should therefore also be taken into 
account^. However, the above discussion on the leading 
effect of t' might leave one wondering whether a pertur- 
bative spin-only Hamiltonian, starting from a half-filled 
one-band Hubbard model can provide a quantitative mi- 
croscopic description of La2Cu04 at all. It is this very 
question that has motivated us in pursuing the work re- 
ported here, and it is why we have derived the spin-only 
Hamiltonian from & t — t'—t" — U Hubbard model, in- 
cluding all ring exchanges and hopping processes to order 
{i/UY- This procedure generates a large number of fur- 
ther neighbor and ring exchange paths, a small number 
of which make significant contributions. The main cor- 
rections come from terms of order {t'Y /U and [t")'^ /U 
but we also find that ring exchanges of order t't"t^ /U^ 
are also significant. The latter, in particular, determine 
the sign of t' /t" . Of course, the resulting theory has now 
two more free parameters (i.e. t' and t") than the orig- 
inal model^ to fit the inelastic neutron scattering data 
and the reader may therefore not be surprised that we 
achieve a better fit than in Ref. However, as we 

show in Section HI, the ensemble of parameters giving 
the best unbiased fit to the data is in good agreement 

with the values suggested from other sources for various 
cupratesia^iSi^i^edo^ 

The rest of the paper is organized as follows. In Sec- 
tion II we present the method to obtain a consistent spin 
representation of the original Hubbard model up to third 
neighbor hops. As an application, in Section III, we in- 
vestigate the consequences of applying the method to a 
general t — t' — t" — U model by fitting the magnon dis- 
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persion data of Ref. [a] for La2Cu04. Section III also 
discusses the procedure to calculate the staggered mag- 
netization operator, and shows that our expansion of the 
t — t' — t" — U model gives a value of the sublattice magne- 
tization in good agreement with experiment^. We con- 
clude the paper in Section IV. Wc have included a num- 
ber of appendices to assist the reader with a few tech- 
nical issues. Appendix A gives the various terms that 
contribute to the effective spin Hamiltonian with arbi- 
trary t, t' and t" up to order In Appendix B we 
show that the renormalization factor for the magnetiza- 
tion, coming from charge fluctuations^ is the same in 
our approach as that found in an extension of the mean- 
field Hartree-Fock method of Ref. [5§| . Appendix C gives 
the k dependence of the various terms coming in the 1/5' 
spin-wave calculation. Appendix D discusses the spin- 
wave velocity renormalization factor, Zc(k). Appendix 
E comments on the results of various constrained fits to 
the spin wave energies. 



II. EFFECTIVE SPIN HAMILTONIAN FOR 

THE t-t' -t" -U HUBBARD MODEL 

A derivation of the spin-only effective theory starting 
from the one band t — U Hubbard model has been given 
by several authorsilii^i^i^^. In this section we investi- 
gate the effects of including new parameters t' and t", 
for direct hops between second and third nearest neigh- 
bors, respectively, on the spin-only effective theory. We 
first recall the key steps in the unitary transformation 
method for the t — U Hubbard model. Wc then apply 
the method to the t — t' — t" — U model and obtain the 
modified spin-only Hamiltonian generated by all possible 
virtual electronic paths up to order l/U^, giving terms of 
order t^U^ {t'f/U^, (^")VC^^ t^{t')^/U^ t^t")^/U^ 
t^t't"/U\ and {t't"f/U^. 



Derivation of the spin-only effective 
Hamiltonian 



We begin with a brief review of the derivation of the 
spin-only Hamiltonian of the one band nearest-neighbor 
Hubbard Hamiltonian, Hn: 



Ha = T + V 



(4) 
(5) 



The first term is the kinetic energy term that destroys an 
electron of spin a at nearest neighbor site j and creates 
one on the nearest-neighbor site i. The second term is the 
on-site Coulomb energy U for two electrons with opposite 
spin to be on the same site i and where rii^a = cj ^Ci,^ is 
the occupation number operator at site i. 

As introduced by Harris et al^ and developed further 
by MacDonald et alr^, the transformation relics on the 



separation of the kinetic part T into three terms that 
respectively increase by 1 (Ti), keep constant(To) or de- 
crease by one (T-i) the number of doubly occupied sites. 
We write : 



T 



-t E , = Ti + To + T_ 
To = -tYl ^i.ffcLcj,<T/ij,ff 



T-i = -t^Kscl^Cj 



(6) 
(7) 



(8) 
(9) 



where a stands for up if a is down and for down if a is 
up. This separation comes from multiplying the kinetic 
term on the right by 71^.5 4- /li^g- = 1 and multiplying on 
the left by n^^g- -I- hj^g = 1- 

Applying the unitary transformation e*'^ to H^i leads 
to an effective Hamiltonian, Hs, where double occupancy 
is eliminated perturbatively, order by order in S. Using 
the relation: 

Hi^) = e^^i/HC-^ = ^^^[^5^_^[^5,[^5,i/H]] , 



1! 



2! 



and, defining as in Ref. [IE 



(10) 



T' -^(to1,TO2, ■ • ■ ,TOfe) = T [m] ^ Trn^Trn^ ■ ■ - T^n^, 

(11) 

S is solved for, order by order. Hence, starting from a low 
energy vacuum of singly occupied sites, Hg contains no 
terms that create or annihilate doubly occupied sites up 
to the order for which S has been determined. This leads 
to an expression for Hg to order t{t/U)^ (see Refs. [l^, 
Hi): 



u 



+ i^T(3)(-i,o,l) 



(12) 



+ ^ (tW(-i,i,-i,i) - rW(-i,o,o,i) 
- ir(4'(-i, -1,1,1; 

where the associated generator for the unitary transfor- 
mation is: 

^S^'^ = ^{Ti - r_i) + ^{[Ti,To] - [To,T„i]) 

+ ^ (- [To, [To, Ti]] - [To, [Ti, To]] - [Ti, [Ti, To]] 
- i[T_i,[ro,T_i]] + ^[Ti,[Ti,T_i]] 



[T_i,[ri,T_i]] 



(13) 
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B. Derivation of the spin Hamiltonian. 



Inclusion of t' and t" 



The effective Hamiltonian Hb^'' in Eq. p2|) is still de- 
fined in terms of fermion operators entering the T'*^^ op- 
erators. When focusing on the magnetic properties of the 
half-filled Hubbard model, it is convenient to recast the 
effective Hamiltonian Hs in Eq. p2p in a spin-only nota- 
tion. For this, we use a mappingi^ between the subspace 
of the Hubbard model with singly occupied sites and the 
Hilbert space of a spin S = 1/2 system. The mapping is: 



Hubbard Space 



Spin — Space 



site i 

site i 



(14) 



(k) 

The spin Hamiltonian Hg acting on this space is derived 
from the Hamiltonian acting on the occupation number 
subspace (see Ref. [l^ for more details) . Once it is writ- 
ten in the spin 1/2 basis it can be transformed into an 
explicitly SU(2) invariant form using: 



1 

2^ 



E 




..Wi/W) 
(15) 



where aaj is the Pauli matrix for site p, with electron in 
spin state ap. From this one finds the spin-only Hamilto- 
nian evaluated to third order in t/U for the t~U Hubbard 
modeli^i^^: 



IT 



<<i,j>> 



4*4 



<i,J»> 



(16) 



Si) 



<i,j,k,l> 

+(S, ■ti){sl -Sj) ~ (s, ■ si) (S^ ■ ^)} 

where < i.j >, << i.j >>, and <<< i,j >>> de- 
note sums over first, second and third nearest neighbors, 
respectively and where < i,j,k,l > identifies a square 
plaquette with sites i,j,k,l defining the four corners of 
an elementary square plaquette. 



(k) 

Having reviewed the procedure for deriving Hs for 
nearest-neighbor hopping t only, we now discuss how to 
include second and third nearest neighbor hopping t' and 
t" . We start with an extended Hubbard Hamiltonian: 



Hh = T + T' + T" + V 

ij3:cr i 



(17) 
(18) 



where t' is the hopping constant to the 2"** nearest neigh- 
bor and t" to the 3'"'' nearest neighbor (see Fig. [T|); Jq, 
is the a nearest neighbor of i. As above, we define the 
operators Tm, taking into account t' and t": 



i.j:a 

To = \ ^{-U,j)(Kac{^Cj,crhj,s 



+ 



T-1 = \ ^{-U,j)Ksc\y 



(19) 

(20) 
(21) 



where ii.j — i \i i and j are nearest neighbors, t^.j = t' 



if i and j are second nearest neighbors, and ti 



i" if I 



and J are third nearest neighbors and otherwise. The 
commutation relations in pop do not change and the ex- 
pression for the unitary transformation in (jl3p remains 
formally the same. On the other hand, the resulting spin 
Hamiltonian gets drastically modified by the inclusion 
of t' or t" . In particular, many new 'plaquette or ring ex- 
change terms are generated. The complete expression for 
the spin Hamiltonian is reported in Appendix ([A]) . As 
in previous wor k^^i'^^ , and even more so here, given the 
much increased complexity of the Hamiltonian, a com- 
puter program was used to collect together all the terms 
in the projection of Eq. [15] which ultimately lead to a 
globally SU(2) invariant effective spin-only Hamiltonian 

i/i")(t,t',t",C/). 



III. PHYSICAL RESULTS 
A. Spin wave dispersion 

1. Spin wave calculation 

In order to obtain the magnetic excitation spectra pre- 
sented in Fig. [3] and Fig. [3] we perform a \/S spin wave 
calculation. The spin operators arc written in terms of 
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boson operators through a Holstcin-Primakof&^ 1/S 
expansion for a bipartite Ncel ordered square lattice: 



Sublattice a 



Sublattice b 

=- s + b]b, 



S+ = a/25 - 



a ai ai 



S- = ^2S-h]b,h, 



(22) 



Si 



In reciprocal space, one obtains to order S the follow- 
ing general expression of the spin Hamiltonian ijf^'' ~ 
-^s 0*^ + -^s quad' where H^^ is the classical ground state 
energy and -fff tiuad is given by 



rr(4) 

s,quad / j 



V 

k 



{a\xi\^ + ^k^k) + Sk (^flkfok + akfok 



(23) 

Note that here and elsewhere in the paper, the sums in 
reciprocal space are over the magnetic Brillouin zone, 
which for a square lattice with N spins contains N jl 

sites, quad '^^'^ diagonalizcd through a Bogoliubov 
transformation 



givmg 



flk = Uk^k + Wk/3k 
6k = "k/3k + Wk^k 



Suad = E^k(4"k + /3i/3k + l) 



(24) 



(25) 



and where 



Wk 



Sk^ 



2ek(Ak - ek) 
Ml - Ck 



2ek 



(26) 



The magnons for a given wave vector k are thus two fold 
degenerate with eigenfrequencies^: 



ek = \ A 



Bl 



(27) 



In this expression, the functions and i?k take into 
account all the bilinear and ring exchange interactions 
entering i/i^^ and given in Appendix |X1 Detailed ex- 
pressions for Ak and i3k can be found in Appendix [Cl 

The staggered magnetization operator, is, in ref- 
erence to sublattice a, defined conventionally for a spin- 
only Hamiltonian as 



i=N/2 



Ml 



iN/2) ^ 



(28) 



Rewriting the previous equation in k space and introduc- 
ing the Bogoliubov transformation above, we arrive at a 



standard expression^: 
2 

N 



^^^^ ^ - ^ E "k^k ^S-^ E("k'al.«k + «k'/3k/3,^ 



off — diagonal terms) 



k 

(29) 



At zero temperature, all (akak) and (/3k/5k) arc zero, 
and one has for the ground state expectation value of 
Ml (Aft): 



fk 



* + 2 iv4^2ek 

k 



(30) 



2. Magnon self energy: renormahzation factor Zc (k) 

The 1/S correction to Eq. ([50)) vanishes for the Heisen- 
berg antiferromagnet, with the next corrections appear- 
ing at order 1/5*^ only^. This explains why spin wave 
calculations give a very good estimate of zero point quan- 
tum fluctuations in this model. On the other hand, the 
terms to quartic order in magnon operators in give 
a contribution to the diagonal part of the magnon self 
energy^. This gives rise to a 1/5* correction to the spin- 
wave energies ek in Eq- (|27p . This renormahzation of 
the magnon energy scale is important for the quantita- 
tive comparison between our calculation and the exper- 
imental results of Coldea et al^ and for the subsequent 
determination of the parameters {t,t' ,t" ,U} and hence 
we evaluate it for our effective spin Hamiltonian. 

Expanding the spin Hamiltonian to fourth order in the 
boson operators of the Holstein-Primakoff transforma- 
tion, one finds: 



(4) 

s,quad 



(4) 

s, quart 



(31) 



quart now treated as a first order perturbation to 



the quadratic spin-wave Hamiltonian 



(4) 

s,quad" 



To this 



order the magnon energy is shifted but the eigenvectors 
are unchanged and magnon-magnon interactions are not 
generated. The shift is found by diagonalizing the set of 
2x2 matrices with elements (0|7k-ffi quart 7k where 
(7k, 7k) £ {o^k, /?k}- We hence find a correction Se^ to 
the magnon energy, and a renormalized magnon energy 
e(k): 

ek — > ek = ek + Se^ = ek(l + Ak) = ek^c(k), (32) 

where we refer to Ak as the magnon energy correction 
and Zc(k) the magnon energy renormahzation factor. 
For the Heisenberg model, the product e^Ak corresponds 
to the leading contribution to the magnon self energy cal- 
culated in Ref. To order t'^/U and with t' and t" 
set to zero, Ak is uniform over the whole Brillouin Zone 
and we find: 



Zc(k) 2± 1.1579, 



(33) 
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as obtained in Ref. [43|. Details of this calculation are 
given in Appendix Id 

When non-zero values of t' and t" are included, Ak is 
no longer constant over the Brillouin zone, but the two- 
fold degeneracy for each wave vector k remains. The 
renormalizing factors 2^c(k) are calculated numerically. 
For a given set of t' and t" values, a finite size scaling 
analysis was used to extrapolate the sums over k' in Ap- 
pendix |D] over the Brillouin zone to the thermodynamic 
limit. In Fig. [2] we show the dispersion of 2'c(k) obtained 
for the parameter set {t, t' , t", U) (see Eq. giving the 
best fit to the magnon dispersion datai^ from Coldea et 
al., as described in the next section. For these values the 
mean value of Ak over the Brillouin zone is found to be 



(Z,(k)) ~ 1.219. 



(34) 



have proportionate consequences for the values of U and 
t deduced by comparisons with experiment. 



3. Experimental data and previous results 

The magnon energy spectrum for La2Cu04 obtained 
through inelastic neutron scattering^ is shown in Fig. [3l 
The data follow a trajectory through the Brillouin zone 
covering, in particular, the line along the magnetic 
zone boundary with wave vector in the interval qez = 
[(tt, 0), (7r/2, 7r/2)], where the previously discussed down- 
ward dispersion is manifest. Also included in this figure 
is the fit of Ref. @ using the Hamiltonian p6p . The pa- 
rameters used there to make the fit at a temperature of 
10 K aro^ 



We also calculate the dispersion of Ak over the Brillouin 
zone for this set of parameters, finding: 



a(Ak) 
(Ak)k 



0.33 X IQ- 



(35) 




0.135 ±0.03 
0.30 ±0.02eV 
2.3 ±0.4eV 



(36) 



where a is the standard deviation of the distribution of 
values. Details of the calculation can be found in Ap- 
pendix |dJ 



-1 
1 




t/u = 


0.1 


±0.05 , 


t'/t = 


-0.35 


±0.05 , 


t"/t = 


0.22 


±0.05 , 


t = 


0.35 


eV 


u = 


3.5 


eV 



FIG. 2: (Color online) k dependence of Ak. Data is shown 
over the structural Brillouin zone of the square lattice. As 
can be seen, the dispersion has the symmetry of the zone 
corresponding to the antiferromagnetic order. 

From our results here one can see that the magnon 
energy renormalization Zc(k), calculated to order l/S*, 
is indeed an important element of the quantitative com- 
parison between our microscopic theory and the experi- 
mental results of Coldea et al^. The values obtained for 
Zc(k) change the energy scale of the magnon dispersion 
by around 20% of what would have been obtained if the a 
bare 2c(k) = 1 had been used. The Zc factors therefore 



While the above fit of Coldea et al^ is clearly good, 
the numerical values should be compared with those from 
experiments and band calculations on other cuprate ma- 
terials. For example, values found by fitting ARPES re- 
sults for the half-filled compound Sr2Cu02Cl2 containing 
Cu02 planes are^^i^i^ 



(37) 



One can see that the energy scale, set by U and the ra- 
222tio t/U in Eq. ((37)) are quite different from that found 
by Coldea et al^. Also, while the effects of t' and 
t" are neglected in Ref. [6|, the ARPES results unam- 
biguously show that t' and i" have non-negligible val- 
ues within the Cu02 plane o^^'^° , although band struc- 
^■'^ure calculationa2ii25i4i on LaCu204 and ARPES exper- 
iments on the lightly doped syste m^^i^°i^-'^i^^ do suggest 
somewhat smaller values of t' than those in Eq. ([571) foi' 
Sr2Cu02Cl2. In other words, based on various experi- 
ments on a number of cuprates, one would have "antici- 
pated" a t/U value for La2Cu04 somewhat smaller than 
the value in Eq. ((551) from Ref. @, although the experi- 
mental uncertainty allows for some overlap between the 
two. Most importantly, from a theoretical point of view, 
t/U — 0.135 corresponds to U — 7.4i, which is smaller 
than the tight binding band width of 8t. This value may 
therefore be a little small^, to insure that the mate- 
rial remains within the Mott insulating phase. The error 
bars do allow values up to [/ = 9.5 {t/U = 0.105), which 
would push the model further into the insulating regime. 
Nevertheless, it is quite possible that the analysis of the 
data in Ref. (d) using the model Hs"^^ in Eq. (fTBl) leads 
to an underestimation of U . 
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(tt.O) (7^/2,^/2) (O.H) 



FIG. 3: (Color online) Magnon energy in La2Cu04, as a func- 
tion of wave vector k, across the Brillouin zone at a tempera- 
ture of 10 K. The experimental data (red squares) and fit (full 
line) are from reference^. The {t,t' ,t" ,17} fitting parameters 
are given in (|36p . 




(7r/2,7r/2) (ir, ir; 



(jr,0) (7r/2,ir/2) ((l.{>) 



FIG. 4: (Color online) Magnon energy in La2Cu04, as a func- 
tion of wave vector k, across the Brillouin zone at a temper- 
ature of 10 K. The experimental data (red squares) are from 
reference^. The fit is made using the parameters given by 

m\i. 



From this discussion it is clear that a more detailed 
analysis of the experimental magnon dispersion that ac- 
counts for t' and t" is needed. 



4- Fitting procedure and fit to the experimental data 

Our theory, H^'^^t, t' , t" , U) in Appendix A, now con- 
tains four independent parameters, t,t',t",U which we 
fit to the ensemble of experimental data points from 
Ref. @. As we show in Appendix [Dl Zc(k) is a func- 
tion of t/U, t' /t, t" /t and k. It thus should be calculated 
over the ensemble of data points for each iteration of the 
fitting algorithm. As the expression for Zc(k) contains 
inner sums (X^k') over the Brillouin zone, the fitting pro- 
cedure is time consuming. A finite size scaling analysis 
of the convergence of these sums towards the thermo- 
dynamic limit is discussed in Appendix |D] The fitting 
procedure is as follows. We choose a first point at wave 
vector k with experimental magnon energy ecxp(k) and 
minimize the quantity (5(k) ~ |e(k) — ecxp(k)| with re- 
spect to t,t',t",U. From this first fit we extract a fac- 



tor it^/U that fixes the energy scale for the ensemble 
of points and allows us to write e(k) as {At^ /U) times a 
function oit/U, t' /t, t" jt and k. This energy scaling step 
can be achieved for many values of t/C/, t' jt and t" jt and 
the constraints on these variables are not very high at 
this stage. The best fit parameter set is now found by 
minimizing the ensemble of variables {(5(k)(t//4t^)} with 
respect to t/U,t' /U,t" /U , using least squares fit. Once 
the final value of t/U is established, the value of t, and 
hence C/, t' and t" , are deduced. The values found are: 



(38) 



The error bars, obtained through the calculation of the 
least-squares function x^, represent a deviation of less 
than 5% from the best fit values. 

If the dispersion in Zc(k) is ignored, it only needs to 
be calculated once for each iteration of the ensemble of 
points, thus considerably speeding up the procedure. Ig- 
noring the dispersion and taking Zc(k) = (Zc(k)) = 
constant, the best fit parameter set for this case gives 
the values: 



t/U = 


0.126 


±0.03 , 


t'/t = 


-0.327 


±0.05 , 


t"/t = 


0.153 


±0.05 , 


t = 


0.422 


eV 


u = 


3.34 


eV 



t/U = 


0.121 


±0.03 


t'/t = 


-0.313 


±0.05 


t"/t = 


0.167 


±0.05 


t = 


0.430 


eV 


u =^ 


3.75 


eV 



(39) 



Comparing the two sets, taking the size of the error 
bars into account and considering the fact that Zc(k) 
represents only the leading correction to the magnon en- 
ergies, we conclude that the dispersion of Ak does not 
have a significant effect on the determination of the pa- 
rameters. 

The best fit to the experimental data using param- 
eter set (pS)) is shown in Fig. Comparing with 
Fig. ([3]), one observes that introducing t' and t" leads 
to a clear improvement in the quality of the fit through 
the experimental points. One might have expected this, 
given the increased number of free parameters. However, 
comparing the set (|38p with ([37|l . one can also see that 
the parameters correspond more realistically with those 
thought to describe other cuprates such as, for example, 
SrCu02Cl2. This is particularly the case for the energy 
scales t and U which are found to be larger than those in 
Eq. The ratios t'/t and t"/t are also in fair agree- 
ment with values for other cuprates and, perhaps most 
interesting, the fitting procedure even finds the correct 
signs for these ratios. The ratio t/U is slightly reduced 
compared with (|36p and more in line with a value that one 
would expect for a system in the Mott insulating regime, 
although the change is less spectacular than in the case 
of the absolute energy scales t and U, with the differ- 
ence in the t/U value remaining within the experimental 
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error bars. However, considering the overall comparison 
between calculations and the experiments of Ref. [6j, it 
would appear that adding all terms in t, t' and t" up to 
order does allow a quantitative description of the 

magnon dispersion in La2Cu04, which is an improvement 
over the equivalent procedure which only includes the t 
and U parameters. 



5. Comments on the influence of t' and t" 

Referring to Appendix \^ one can see that the inclu- 
sion of t' and t" and all four hop processes involving t, t' 
and t" introduces a large number of new ring exchange 
terms. From this observation one might ask which of 
the many terms make the most difference and could we 
have got away with a less complex effective spin Hamil- 
tonian? In this section we address these questions and 
justify our choice to do a complete calculation of all four 
hop processes. 

The first question is whether or not we can limit our- 
selves to just one extra parameter t' or t" . This point 
has already been discussed in the Introduction (Section 
II Bl) . We argued there, (Eq. ([3])), that t' increases the 
strength of the antifcrromagnetic second neighbour ex- 
change, J2 t'^/U'^ + 4(t')^/C/. This in turn modi- 
fies the effective second neighbour coupling in the spin 
wave analysis, — > 4(<')^/C/ — IQf^ /U^ making it less 
strongly ferromagnetic, or even antifcrromagnetic if t' is 
sufficiently large. As the negative magnon dispersion 
along the zone boundary^ is a consequence of a ferro- 
magnetic (i.e. negative) Jf'^, one expects that fitting the 
calculated magnon dispersion to the experimental data 
will require a larger value oi t/U if t' ^ but t" is ig- 
nored. This is indeed the case; introducing t' to order 
{t'Y /U and setting t" = 0, as well as neglecting all terms 
of order {t'tf /U^ and (^')VC^^ the following best fit pa- 
rameters are found: 

r t/U - 0.166 ±0.016 , \ 
\ t'/t = -0.24 ±0.01 , J ■ 

The results are very close to those found from a random 
phase approximation^^i^ and also from quantum Monte 
Carlo calculations^. It is clear that the addition of t' 
only has taken the parameter set in the wrong direction: 
t/U is now further from most estimates than the case 
when t' = 0. Also, with U = 6t the system could be very 
near the metallic phased, and not deep in the Mott insu- 
lating phase, as is the case for La2Cu04. As we discussed 
in the introduction, for U = 6t the t' = model is right 
at the metal-insulator transition and the presence of a 
finite t' increases even further the value of U necessary 
to position the system in the insulating phas o^^i^^ . 

The leading effect of third-neighbor hopping t" is also 
to increase an exchange interaction, this time J3, through 
the relation : 



•^3 ^3 =4^+4i^. (41) 

However, in this case antifcrromagnetic (positive) ex- 
change also contributes, as does the ring exchange term, 
to a negative (downward) dispersion along the magnetic 
zone boundary. Including just t" , one would therefore 
expect the best fit parameters for the data to include a 
particularly small ratio of t/U , which is indeed the case. 
Hence we see that, as far as this characteristic dispersion 
is concerned, t' and t" play opposing roles. It is therefore 
clear that one needs to include both hopping constants 
to get a good fit to the experimental data as we have 
illustrated in Fig. 2] with parameters (Eq. (P7| ) similar 
to those in other cuprates^^i^liii. 

Given that individually, terms of order t^t'^/U^ are 
only one tenth of the magnitude of terms of order t''^ /U 
or of the terms of order t^ /U^ , one might expect the 
good fit reported in Sectionf [38)) to be also obtained when 
taking into account terms of order {t'Y /U and {t")'^ /U 
only, and ignoring the host of more complex four hop 
terms generated by the further neighbor hopping. Not 
so! If one neglects them and only considers t' and t" to 
order A{t'Y /U and A{t"Y /U), one obtains the following 
set of best fit parameters: 



t/U 


= 0.14 


±0.03 , ^ 


t'/t 


= 


±0.0, 


t"/t 


= 


±0.0 , ► 


t 


= 0.303 


eV, 


u 


= 2.14 


cV 



Within numerical error this is the same data set as 
Coldea et al.{ [SB]), with t'/t = t" /t = 0. The rather sur- 
prising and interesting conclusion is therefore that, on 
their own, the two-hop further neighbor terms make no 
net contribution. Rather, it is the four-hop processes that 
make the difference. It seems that their small value is 
compensated for by their multiplicity; that is, by the fact 
that a given site i appears in a large number of diagrams 
of order 1/U^. Hence, their contribution is increased by 
an order of magnitude and in this intermediate coupling 
regime, the global effect of the near-neighbor four hop 
terms, the majority of which are ring exchanges, is as im- 
portant as that for the two hop further neighbor processes 
(Note that as t' and t" arc both zero in Eq.lgll), Zc(k) 
is once again given by (|33p over the whole zone). An 
alternative procedure would be to impose some partial 
constraints on the parameters {t,t' ,t" ,U), using values 
obtained from, say, band structure calculations. Results 
of such fits are briefly discussed in Appendix [El 

Finally, we remark that in our fitting procedure, the 
symmetry is broken between data sets with positive and 
negative values for the ratio t'/t". The best fit occurs for 
a negative ratio t'/t" (see Eq. (p7|) '). Recall that a canon- 
ical particle-hole transformation at fixed chemical poten- 
tial changes the sign of all hopping integrals and leaves 
the filling invariant, since 2 — n = n when n = 1— An 
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additional change of phase of the creation-annihilation 
operators on one of the sublattices by tt restores the sign 
of nearest-neighbor hopping t, but not that of t' and t". 
This shows that t'/t and t" /t can have arbitrary sign but 
that the sign oit'/t" is physically relevant. The symme- 
try breaking between positive and negative signs for t'/t" 
can only be accessed by including four hop processes. As 
one can see from Appendix |X1 most of the spin interac- 
tion terms are even in powers of t, t' and t" . However, 
a small number of terms have odd powers of t' and t". 
These are the terms that determine the sign of the hop- 
ping ratios. This can be seen in Fig. [5] where we show 
the parameter from the fitting algorithm as a func- 
tion of t'/t and t" /t (the value oi t/U being set to the 
best fit value). If the terms proportional to t^t't" /U^ are 
r(4) 



deleted from Hs (see Fig. 5(a) I, the x parameter has 



four fold symmetry and there are four points of best fit 
independently of the sign of t'/t and t" /t. On the other 
hand, when these terms are included (see Fig. |5(b)"] ), the 
symmetry is clearly broken. The minima for positive ra- 
tio t'/t" become broader and shallower, while those for 
negative sign arc narrower and deeper. The fitting pro- 
cedure hence favors the two points with a negative sign 
for the ratio t' /t" . To lift the remaining degeneracy we 
have to appeal to band structure calculations^ii^iii or 
ARPES experimentsiai^i2ii22^^ l^o^j^ of ^l^j^.;^ ^y^ 

negative, hence t" /t positive. 



B. Magnetization and spin Hamiltonian 

Having found a set of parameters {t,t' ,t" ,U) that 
suitably describes the experimental spin-wave dispersion 
data, we now explore how this new set and the derivation 
of the effective spin Hamiltonian, Hs^\t,t' ,t" , U) affects 
the value of the zero temperature Neel order parameter. 




(a)x^ with odd powers of t' /t and t"/i 
excluded 




(b)x^ with odd powers of t'/t and t" /t 
included 

FIG. 5: (Color online) Goodness of fit parameter (grey 
scale) as a function t"/t and t'/ t wit h the terms wi th od d 
powers of t'/t and t" /t excluded (5(a) I and included (5(b) I 



1. Magnetization operator 

In the construction of effective theories all operators. 
Oh, defined in the original model must be canonically 
transformed before they can be exploited in a calcula- 
tion within the spin only theory. That is, within the 
effective theory. Oh becomes Og = e^'^One"*'^ and the 
expectation value in the ground state is defined by 

^ h(0|Qh|0)h ^ s(0|Qs|0)s 

h(0|0)h .(0|0)3 ■ ^ 

Here |0)h and |0)s = e*'^|0)H arc the ground state wave 
vectors in the original Hubbard (H) and spin-only (S) 
models. We have recently showii^^ that incorporating 
this transformation has important consequences for the 
ground state magnetization as one moves into the inter- 
mediate coupling regime. On application of this proce- 
dure on the magnetization operator we find^ new quan- 
tum fluctuations arising from the charge delocalization 



over closed virtual loops of electronic hops. These spin 
independent fluctuations allow us to reconcile an appar- 
ent paradox concerning the behavior of the ground state 
magnetization, as a function of t/U. As discussed above, 
the leading effect of including processes to order t{t/U)^ 
is to introduce effective second and third neighbor spin 
interactions (see Eq. ((!])). Within the lowest 1/5' order 
spin wave approximation, the resulting effective ferro- 
magnetic second neighbor interaction reduces the trans- 
verse quantum spin fluctuations and stabilizes the Neel 
order of the spin model. This would naively seem to 
imply that the staggered magnetization should be an 
increasing function of t/U as the system departs from 
the Heisenberg limit, t/U 0; a result which is diffi- 
cult to justify on physical grounds. Indeed, one might 
expect that as t/U is increased, the enhanced electron 
mobility would lead to a progressive return to the non- 
magnetic metallic state, with the staggered moment de- 
creasing as t/U increases, as found in Ref. [5^, rather 
than increasing. In Ref. |33| we showed that this is indeed 
the case: the new quantum fluctuations, or "charge flue- 
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tuation" channel, arising from virtual doubly occupied 
states, counter the effect of generating a ferromagnetic 
second neighbor interaction. Consequently, the ground 
state magnetization for the Hubbard model, when calcu- 
lated using the spin-only description, is indeed a decreas- 
ing function of t/U. The main steps of the calculation are 
reviewed below and then extended to include the hopping 
parameters t' and t" . 

Within the Hubbard Hilbert space, the staggered mag- 
netization operator is defined by: 



numbers of hops, such as the second term in the right- 
hand side of Eq. ([5T|) . are zero. Writing the expression 
in terms of S* = 1/2 spin operators^ we find: 



Ml = 



<»,i> 



(52) 



H - 



-iQri 



(44) 



In the case of periodic boundaries, or in the thermody- 
namic limit where boundaries can be neglected, all sites 
become equivalent and the correction to the magnetiza- 
tion operator becomes a multiplicative factor—: 



where: 



(53) 



Q = (tt, 7r) and 5',f = - ^ 4 si^si.s2'^'^s2 • (45) where Ml = ■j^Ylii^l^''^"'" is the naive expression for 



The magnetization in the ground state |0}^ of the Hub- 
bard model is defined as: 



A/ = h(0|a4|0)h. 



(46) 



Working within the effective spin theory, we have to ex- 
press all the operators in the spin language. Mathemati- 
cally this means that the unitary transformation in ([1 
has to be applied to all operators including : 



(47) 



The commutation relations between M^j and the dif- 



ferent operators in ()13p an 



t 



N 



N 



(48) 



TV 



N 



(49) 



1 ^ 



K^cl^c,^,K-,{-irK,a) ■ (50) 



Expressed in terms of these operators, the staggered mag- 
netization operator reads^: 



Ml = M + - [f^-f^i 



1 

2C/2 



2 ,T^iTi-T^iTi 



(51) 

The symmetry associated with half filling ensures that 
contributions to this expression corresponding to odd 



,33 



^ TV 

the magnetization operator in the spin languag( 



2. Role of t' and t" 



Inclusion of the higher order hopping constants t' and 
t" introduces extra paths for charge delocalization and 
we might expect new charge fiuctuation terms to appear. 
The hopping operators Tm in Eqs. (|19l20l2ip . can be 
generalized to include t' and t" and the commutation 
relations in Eqs. (|48l49l50p are redefined: 



[Tm, Mjj] — j^Tm, 



(54) 



(55) 



(56) 



for hops involving and t", respectively. The stag- 
gered magnetization operator (|5H) is thus also general- 
ized. Dropping the first order terms in \/U that do not 
contribute to the expectation value in a singly-occupied 
state, Ml takes the form: 



Ml 



M 



^ { (f_i + f '_i + f"_i) (Ti +Ti+ T[) 



i')} 



. (57) 

However, it turns out that the operators T^j and de- 
fined in Eqs. (|55l56p are identically equal to — . From 
this result, it follows that the first corrections to Eq. ((53|) 
that depend on t' and t" appear beyond the second order 
in the perturbation scheme. Hence we find that to order 
(t^/C/)^ the expression for Ml{t,t' ,t") is unchanged and 
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given by Eq. (|53|) . However, this does not mean that 
introducing t' and t" has no effect on the expectation 
value for the staggered magnetization. Firstly, the in- 
clusion of t' and t" in the magnon dispersion modifies 
the ensemble of parameters {t,t' ,t" ,U) that best fit the 
data. Secondly, the zero point fluctuations are controlled 
by {t/U,t' /t,t" /t), through the k dependence of Uk, Wk, 
and i?k in Section If 11 A II and this, as we discuss 
below, directly affects M^. 

The charge renormalization factor has also been cal- 
culated using a Hartree-Fock mean-field method as in 
Ref. [5^, in which transverse spin fluctuations are ne- 
glected (see Appendix B). This calculation on a square 
lattice gives the same factor, 1 — 8{t/U)^, independently 
of t' and t" . Physically this comes about because, in 
the spin-density wave ground state (i.e. Necl order, 
Q = (7r,7r)), electron hops between two sites that be- 
long to the same Neel ordered sublatticc arc prohibited 
by the Pauli principle. Hence, same-sublattice hops by 
t' and t" do not contribute to charge renormalization to 
order {t')'^/U, {t"f/U. 



3. Consequences 

An experimental estimate of the zero temperature sub- 
lattice ordered moment in La2Cu04 has been determined 
from neutron scattering. Lee et al^ find 

Mexporimcnt = 9 {MI) l^B = 0.55 ± 0.05^*5 , (58) 

where g is the Cu^+ Lande p-factor. The p-factor for 
Cu^+ in layered cuprates such as La2Cu04 has appar- 
ently not been determined experimentally^. However, 
a typical value for a distorted octahedral environment 
is given by Abragam and Bleaney^, g ~ 2.2. We use 
this value to compare theoretical estimates of the or- 
dered moment with the experimental value. A spin wave 
calculation for the nearest-neighbor Heisenberg antifer- 
romagnet gives^i^i^ (Mgt)"" ~ 0.304, which is in good 
agreement with the value obtained from quantum Monte 
Carlo. This gives for the moment, IVF"" = g(Mj'^}"'^, 
M^^ K, 0.67/iB, well in excess of the above 0.55/^b value 
determined experimentally^. Using the parameter set 
of Coldea et al^ but neglecting charge renormalization 
leads to an increased expectation value: — 0.32, 

taking the moment to 0.70/1^, hence further from the 
experimental estimate of Ref. [4J| . Including the charge 
renormalization factor— 1 — 8t^/J7^ of Eq. (|53p gives 
{Ml) ~ 0.27 and M^*) ~ 0.59/iB in better agreement 
with experiment. Including both the charge renormal- 
ization and the effects of t' and t" we find 

(Mt)^"'*"' ~ 0.235, (59) 

from which we obtain 

(Mt)("'*"Ug(Aft)("'*")^ 0.52ms , (60) 



which is in even better agreement with experiment. 
Given the uncertainty in the experimental measurement 
value in Eq. ((571) ^^d in the theoretical value that should 
be taken for g, it is difficult to make in-depth comment on 
the difference between (Af)** * and (A/)*. However, the 
results do serve to illustrate the importance of the charge 
renormalization at this level of approximation. Here, 
it is an essential element in the quantitative agreement 
with experiment, and it is only through the inclusion of 
spin independent quantum corrections (e.g. virtual dou- 
ble site-occupancy) that quantitative agreement can be 
achieved. 



IV. DISCUSSION 

It is known from ARPES measurements and band 
structure calculations that t' and t" are not negligible 
compared to i in a number of copper oxide materials. We 
have thus introduced second and third neighbor hopping 
parameters, t' and t" , into the Hubbard Hamiltonian and 
included all closed four hop virtual electronic pathways 
into the canonical transformation to derive an effective 
spin-only Hamiltonian for the case of a half filled band. 
This Hamiltonian contains many ring exchange terms. 
These kinds of terms have been the subject of much dis- 
cussion Iately^i24i^i2&i^i2^. Through this calculation we 
are able to test the capacity of the one band Hubbard 
model to describe the magnetic properties of the anti- 
fcrromagnetic parent high-temperature superconductor 
La2Cu04. We find in general good quantitative agree- 
ment between the predictions of the t — t' — t" — U model 
and experimental measurements of the magnetic proper- 
ties of the compound^. 

An exact solution of the spin-only Hamiltonian we have 
derived is clearly far out of reach. In fact, because of the 
frustration and ensuing sign problem introduced by the 
various frustrating bilinear spin-spin couplings and the 
many ring exchange terms, even a quantum Monte Carlo 
attack on the low temperature properties of this model 
seems difficult to envisage in the near future. However, 
it has long been known that a spin-wave analysis up to 
leading order in a \/ S expansion accurately reproduces 
the zero temperature staggered moment in the Heisen- 
berg model on a square lattice^. This is because of the 
cancellation of the magnon-magnon interaction terms to 
second order in l/S" for this particular lattice^. This sug- 
gested that a similar spin wave analysis can be a good 
starting point to describe the observed magnetic excita- 
tion spectrum in La2Cu04 and this is indeed what we 
have found in this paper. 

In our spin wave calculations, we calculated the mag- 
netization to leading order in 1/S and included the first 
order correction to the classical spin wave frequencies. 
The magnon energies are thus renormalized by a spin- 
wave renormalization factor Zc(k). Since the classical 
spin- wave frequencies are already of order we only 

retained terms of order /U, {t'Y /U and {t"Y /U in cal- 
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culating Zd'k). Hence, to this order, it does not depend 
on the ring exchange terms which are of order In 
the end, Zc(k) raises the energy scale by about 20% over 
the whole of the Brillouin zone, allowing good agreement 
between experiment and theory without the introduction 
of any arbitrary scale factors. 

To make contact with experiments on La2Cu04 the 
spin wave calculation was set up starting from the clas- 
sical antiferromagnetic Neel ground state. Within this 
framework, the general effect of the ring exchange is to 
reduce the transverse quantum spin fluctuations rather 
than increase them (see however the discussion of charge 
renormalization in the following paragraph). With the 
parameters coming out of the calculation, we therefore 
seem far from entering an exotic quantum phase driven 
from large quantum spin fluctuations about the Neel or- 
dered state. Yet, the ring exchanges arc found to play 
a very important role. In particular, it seems that the 
large numbers of such terms and their large dimension- 
less prefactors (see Appendix A) outweighs their small 
0{1/U^) amplitude, although a more extensive study is 
required here to understand this point in detail. Further, 
ring exchange plays a crucial role in determining the rela- 
tive sign of the ratios t' /t and t" /t. There is experimental 
evidence from ARPES measurements and band structure 
calculations that t' /t is negative while t" /t is positive in a 
number of cuprates. In our calculation there are terms of 
order t^t't" /U^ which break the t" /t' -t" /t' symme- 
try in the effective spin-only Hamiltonian ifi^^ , allowing 
us to deduce the correct sign for the ratio of t' /t" com- 
pared with experiment. In our procedure, we determined 
unconstrained values of U, t, t' and t" by comparing our 
magnon dispersion data with that from experiments on 
La2Cu04^. The best fit is found with one positive {t" /t) 
and one negative {t" /t') ratio, exactly as in ARPES ex- 
periments and band calculations and, for example, our 
best fit data set is in good agreement with those from the 
ARPES experiments on Sr2Cu02Cl2 (compare Eq. pS]) 
with Eg. ((57)) ). 

At the level of approximation at which the canonical 
transformation is performc d^^'^^i'^^ , finite charge mobil- 
ity rcnormalizcs the magnetic moment calculated from 
the spin only Hamiltonian by a factor (1 — 8t'^/U) — . 
We found that t' and t" do not contribute further to this 
factor. This is because direct second and third neigh- 
bor hops take electrons from one site to another on the 
same sublattice, hence virtual double occupancy is ex- 
cluded by the Pauli exclusion principle. The expectation 
value for the zero temperature staggered moment coming 
from the spin only Hamiltonian is thus scaled by the same 
{l-8t^/U) factor as in the case where t' t" = 0^. The 
relevance of this term can be estimated by comparing the- 
oretical and experimental estimates of the total magnetic 
moment. Here, we also find good agreement between the 
estimates of Lee et al in Ref. and our calculated value. 
The charge renormalization is important here, as it scales 
the magnetization obtained by considering solely trans- 
verse spin fluctuations, by an amount clearly in excess 



of the experimental error bars. The comparison with ex- 
periment therefore provides an implicit illustration of its 
importance in the real material. 

Our results also bear on the question of the appropriate 
ratio of the bandwidth ^ it to the interaction strength 
U. First, note that on purely theoretical grounds, when 
U is of the order of the tight-binding band width, one 
expects the Hubbard model to undergo an insulator to 
metal transition, resulting in a breakdown of the per- 
turbation expansion in t/U ^L^, In our calculations, 
since the next order terms in the expansion of t/U is 
smaller than the leading term, there is no evidence of 
such a breakdown. In fact, for the value t/U 2± 0.126 (or 
U 2± 8t) that wc found, we argue that the Hubbard model 
is definitely in the Mott insulator regime, as is necessary 
for the expansion to be valid. Indeed, at ti = t" = 0, one 
can extract from quantum Monte Carlo calculations^ 
that the Mott gap should close around t/U ~ 0.167 (or 
U ~ 6t), a value consistent with recent estimates from 
quantum cluster calculation a^^'^^ . Although the criti- 
cal t/U in general depends on t/^i^ for tf = — 0.3t on 
the square lattice, the critical t/U is still very close to 
the value appropriate for tf = 0^. One can also check 
that for t/U ~ 0.126, the single-particle spectral weight 
displays bands associated with antiferromagnetic excita- 
tions that disperse with J and are distinct from Hub- 
bard bands further away from the Fermi energyi^ When 
the gap is induced purely antiferromagnetic fluctuations 
(Slater mechanism), the distinct Hubbard bands arc ab- 
sentJ^. In addition, for insulating behavior induced by 
antiferromagnetic fluctuations, the potential energy de- 
creases when insulating behavior occurs This is not 
obseved for t/U ~ 0.126^. 

For t/ = t" = 0, the best fit value of t/U ~ 0.14 found 
by Coldea et al^ places La2Cu04 dangerously close to 
the insulator to metal transition discussed above. How- 
ever, with the inclusion of t' and t" and of all ring ex- 
change terms to order 1/U^, we find a ratio t/U ~ 0.126 
which is decreased compared with the initial fit of Coldea 
et ai, which did not include these further neighbor hop- 
pings. This places the ratio t/U for La2Cu04 within the 
Mott insulating phase of the Hubbard model discussed 
in the previous paragraph, by contrast with the result 
found in Ref. HJ. If we had obtained the opposite re- 
sult, namely that the best fit for t/U increases in the 
presence of t' and t" , the whole approach would have be- 
come extremely doubtful as a description of La2Cu04. 
Another approach, starting for example from the three- 
band model^, would have become necessary. Instead, 
the result that we flnd here for t/U gives a consistent 
picture where the parameters of the one-band Hubbard 
model describe an insulator at half-flUing. These parame- 
ters are also in agreement with those describing the doped 
insulator—, including the sign and magnitude of the ratio 
t' /t" . In fact, a non- vanishing t" is necessary to obtain a 
consistent picture since, without it, t/U would increase^ 
compared with its t' = t" ~ value, an undesirable state 
of affairs, as argued above. 



To close, and reiterate, one of the noteworthy results of 
this work is that the presence of the ring exchange terms 
in the effective spin-only theory allow us to determine 
the relative sign of t' and t" . More generally, the theory 
presented herein, restricted to the reduced Hilbert space 
of spin degrees of freedom only, gives results for the elec- 
tronic band parameters that are compatible with ARPES 
experiments and theoretical band calculations for a vari- 
ety of cuprate materials. In addition, the value of f7 ~ 8t 
that we find is consistent with the origin of the insulat- 
ing behavior in parent high-temperature superconductors 
being mostly due to Mott Physics. As a final word, we 
propose that our results provide further evidence that the 
t—t'—t" — U one-band Hubbard model in the intermediate 
coupling regime gives a consistent unified description of 
high-temperature cuprate superconductors and of their 
parent insulating phases. 
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Figure [5] illustrates the sites involved in t' and t" hopping 



i j m 



1 k 

FIG. 6: Label of the different sites involved in the many-spin 
terms induced by t' and t" to order l/U^ 



(sec also Fig. ([IJ). The symbol — X means that the 

corresponding sites participate in the expression of the 
coupling interaction, whereas the symbol — (J — means 
that these sites are transparent in the electronic pro- 
cess. That is the electron 'hops over' the site in question. 
As there is at present considerable interest in diluted 
Mott-Hubbard systems, such as in the La2Cuj^Zni_^04/ 
La2Cu2;Mgi_a:04/iSiZiiSiZ&j we keep track, in the deriva- 
tion of the spin Hamiltonian, of the occupation of the 
sites visited by the electrons, so that ei = 1 if a site i 
is occupied by a spin and e; = if the site is not occu- 
pied. Here we only consider the case of 100% Cu^+ site 
occupation (x = 1). 



APPENDIX A: SPIN HAMILTONIAN WITH t, t' 

AND t" 

The purpose of this appendix is to give the new spin- 
spin interaction terms appearing in the spin-only Hamil- 
tonian including t' and t" up to order l/U"^ in the canon- 
ical transformation. The sole effect of introducing t' and 
t" to order {t'Y /U and {t")"^ /U is to rcnormalize the cou- 
pling constants already present in the spin Hamiltonian 
(Eq. [TB]) generated from the t — U Hubbard model, up to 
order t^/U^. Including t' and t" up to order, \/U^ gen- 
erates many more terms that further renormalize these 
coupling constants. However, a number of terms with 
new "topologies" appear and these play a key role on the 
results discussed in this paper. For reasons of compact- 
ness we introduce the following notation: 



X — yr- 



¥ — ¥ 



>^ — ^ 

1 k 



¥ — ¥ 



>^ — ^ 

1 k 



A — Si - Sk {ejEfc} 

4 Si ■ Sm {CiEjn} 



t'h^ 
C/3 



4 -Jp- {e^£j£k} {S^ ■ Sj + Sj ■ Sk) 



im-^{S,-Su){SySi) {e.e.efeQ} 
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80-^ Ciejekei Pi 



i.j,k,l 



4-^^ Si • Sn {^ii-j^n} 



Si ■ Sn {^i^j^n} 
on -D^-J^^^n 

4-^e»ejemen(-l) P2 



80— i'i"e.e,e™e„ Pi"''"''" 



4— f i" e.e,e™e„ (-1) P^'^'"'™ 



f/3 



80^:^ €i€jemik Pf 



i.j,k,m 







k 



i'^— oJ( k 

i)L-oJ( k 
'Sv^k 



i,j.k,l 



80-^e^eje,efc Pi 



160 ^3 f-if-jeiet {{Si ■ Sk){Sj ■ Si)} 
+80— e^eJelek P'{'''''^ 



t t" 

4 —TTTT- eiejeiek (-1) P2 



i,j.k.l 



C/3 



C/3 



."4 



i.j,k,l 



fff 



2^/2 



i.j,k,l 



APPENDIX B: MEAN-FIELD SOLUTION OF 
SPIN HAMILTONIAN WITH t, t' AND t" 

The purpose of this appendix is to determine the 
sublattice magnetization at zero temperature for the 
t — t' — t" — U Hubbard model using the Hartree-Fock 
(mean-field) method of Ref. [H. From Eqs. pTTT51) . the 
t-t' -t" -U Hubbard model is: 



Hn = T + T' + T" + V 



(Bl) 



(B2) 



4-^ Ciej-e^efc (-1) P^' 



where ji, j2 and js are respectively the first, second and 
third nearest neighbors of i. Fourier transforming this 
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expression leads to: 



M = {n\Si^\n) = ^^^^k^^k ; 



k,(7 
U 

2N 



hence: 



X! X! '^'T,<T'^,3,/3'4',fTcLk'+q,/3'C-k+q,/3Ck,, 



M 



k,k',qcr,(T',/3,,3' 



where: 



ek 



— 2i(cos(fc:c) + cos(fcj^)), 

—2t' {cos{kx + fcy) + cos(fc3; — fcj,)), 

-2t"(cos(2fc^) + cos(2fcy)). 



(B3) 



(B4) 



The sublattice magnetization for nesting wave vector 
Q is defined by: 



M - {n\si^\n), 



(B5) 



where we consider the case where the magnetization is 
polarized along the z direction. \^^) is the spin density 
wave ground state and Q. The spin density operator 5q 
for arbitrary q wave vector is defined by: 



(B6) 



k,Q,/3 



N being the number of sites. 

We define new operators 7^ ^ and 7^ ^ through the 
Bogoliubov transformation: 



7k,a 



MkCk.a + ^^k ^ 0-Q,/3Ck+Q,/3, 
fkCk.a ~ Mk ^ (jj^^^Ck+Q,/?- 



(B7) 



and inject (|B7p into the Hartree-Fock factorization^- of 
Hamiltonian (jB3p . finding: 



E^ = e'^ + e'i± VA^ + ek^, (B8) 
where the density wave gap A is given by: 

A = m 

The Bogoliubov coefficients for this transformation Uk 
and Uk are: 



Wk' = ^ I 1 



ek 



I'k 



fk 



i?k-el 



(BIO) 



With the valence band filled and the conduction band 
empty, since the calculation here is done at zero temper- 
ature, it follows from the above relationships that 



y e . 



(Bll) 



(B12) 



Putting (|B9p into (|B12p gives the following self-consistent 
equation for the sublattice magnetization: 



U N ^ 



'Ck^ 



(B13) 



We remark that equation (|B13[) is independent of ej^ or 
e'^. Considering t/U 1, we can expand the square 
roots of (jB13p in cu/U, after which wc find: 



-E- 



iV^(ek2 + A2)i/2 

k 



-E- 



N^UM 

k 



Further simplifications lead to: 
A/ = 1 4 



^E 



1 - 



ek 



k 



U^APJ U' 
(B14) 

(B15) 



To first order in t/U we obtain 
M 



1 - -E^- 



Since ^E^^ = -^^EJcosCfc.) + 



cos 



(B16) 



i{ky)f, and 



^5](cos(fc,) + cos(fc,))2 = i (B17) 



we 



conclude that: 



M 



'TP- 



(BI8) 



This result is equivalent to that obtained through the 
unitary transformation method in Eq. (|53[) if trans- 
verse spin fiuctuations are neglected, as they are in the 
Hartree-Fock approach. Mathematically, it is the prop- 
erty Ek = — Ek-i-Q that makes nearest-neighbor hopping 
different from the other two hops which both obey, in- 
stead, e'^ = ek_|_Q, ^k = ^k+Q- Note that ek, ej^., e'^ in this 
appendix arc the quasi-particle energies, which should 
not be confused with the spin- wave (magnon) energies in 
Section UTTXTl 

Physically, since t' and t" are hops within the same 
sublattice (parallel spins) of the classical Ncel antiferro- 
magnetic solution, in the strong coupling limit the Pauli 
principle prohibits hopping between these sites. Hence, 
to leading order, t' and t" cannot change double occu- 
pancy or, in other words, produce charge fluctuations 
(or "charge renormalization" ) . 
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APPENDIX C: SPIN WAVE RESULTS space, and collecting all terms, we obtain, after some te- 

dious but straightforward algebra, the following expres- 
By introducing the transformations of Eq. ^ in sions for the functions At and Br: 
-^s quad(^' ^) (cquation (HJ)), passing into Fourier 
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U 

W) 



Ak = 4Ji5 + J2S 



4 cos{kx) cos{ky) — 4 



JaS 



2cos(2fc2:) + 2cos(2A;y) -4 



4 + 4cos(fc^) cos(fcy) 

— 4 + 2cos(fci, + fcj,) + 2 003(^3; — ky) 
16 — 16cos(fca;) cos(/cy) + 4cos(2/ca;) + 4cos(2fey) 
cos(A:2;) cos(fcy) — 1) + 4(cos(2A:^) + cos(2fcj,)) + 16 + 8(cos(fca;) cos(fcj,) — 1) 



U J \t 
2 



+ 125 
+ 2053 



2 /^„x 2 



t \ ^ ft"^ ^ 



-2053 ( _ 



- 16 + 8cos(2fc^) + 8cos(2fcj;) 



32 + 8(cos(2A:^) + cos(2fcy) - 2) 



t'\ ft 



t \ t 



t'\ ft" 



t 

U J \t J \ t 

2 



s 



32 — 8{cos{2kx) + cos{kx — ky) + cos(fci: + ky) + cos(2fey)) 
32 + 16(cos(fc^) cos(fcy)) + 8(cos(2fc2;) + cos(2fcj^)) 

16 - (16(cos(A:^) cos{ky) ~ 1)) - 8(cos(2A:j;) + cos(2fcj^) - 2)) 
+ 4(4cos(fci^) cos(fc2^) -4) 



U 



2053 



t'\ ft 



t \ t 



s 



t ft'V ft" 



cos(fca;) cos(fcj,) + 4cos(2fca;) + 4cos(2fcj,) 

2 

. 1^; \ f I —32 + 16(cos(2fca;) + cos(2fcj,)) + 16 cos(fca;) cos(fcj, 
- 8cos(2fc^ + 2ky) - 8cos(2fc3; - 2ky) 
25 ( (t) (t) -48 + 8(cos(2fc:,) + cos(2fcy)) + 24(cos(fc^) cos(/cj,)) 

+ 4(cos(2fc;r + 2ky) + cos(2fc2; - 2ky)) 



S 



t \^ ft'V ft"^ ^ 



U 



t \ ft"-- ^ 



tV ft'^' 



-s 

4053 



u J \t 
t'V f t^'^ 



32 cos{kx) cos(fey) — 32 

— 16 + 4(cos(3fca; + ky) + cos(3fca; — ky) + cos{kx + 3fey) + cos(fc2: — 3fcy)) 
- 4 + 4(cos(2fc^) + cos(2fcy)) - 2cos(2fc^ + 2fc.y) - 2cos(2fcj, + 2ky) 
{4:COs{2kx) + 4cos(2A:y)) - 12 + 2cos((2A;2; + 2ky) + cos(2fc2: + 2fcy)) 
16(cos(fcj;) cos(fcj,) — 1) 



16 + 16cos(fca;) cos(fc.y) 



(CI) 
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2JiS 



COs(fca.) 



- cos(fcy) 



COS 



\{kx) + cos(fcy) 



1-5 (yjj^ (^"^) ^ cos{2kx + fcj,) + cos(2fcj, + k^) + cos(2A:2; — ky) + cos(2fcy — fca;) 



+80^3 J- _j cos(fc,) + cos(fcy) 



2 / ,/\ 2 



t\ ft' 



2 / 



tY ft' 



4(cOs(fca;) + COS 

2(cos(3fca;) + 



cos(3A:.y)) 
+ 2 cos(fcj; + 2ky) + 2 cos(fcj; — 2ky) 



iiky)) + 4(cos(A:a;) + cos(fcy)) + 4(cos(A;a;) + cos{ky)) 
2 cos(2/c.r + ky 



- 2 cos(2fc^ — kyj 



h20S^ I 



- 8(cOs(fca;) + COS 



iky)) 



- 4(cos(2fc 



fcy) + COS(2fc3, — l\,yj 



— ky) + COS(fc; 

8(cos(fc^) + COS 



- 2ky) + cos{kj. 



- 2ky)) 



- ky) 



- co's,(kx 



- 2ky) + cos(fc: 



- C0S(2fcy + fca;) + 
< \ ^ ft2\ ft^ 



t \t 



i{ky)) + 4(cos(2fca. + ky) + cos(2fc.r 

- 2ky)) 

A{2cos{kx) + 2cos{ky j) - 4(cos(2A:a; + ky) + cos(2fca; - ky) 
cos(2fcy ~ kx)) 
24(cos(fca;) - 



' cos(2fc.y + fcj:) + 



COS(2fc.y 



- cos(A:y)) 



- 4(C0S(2fc3; + fcy) + C0S(2fca; — fcy ) 



^ ft' 



^ cos(2fc. 



fcy) + COS(fca; 



8(cos(fc:r) + 

2fcy) + COs(fca: 



cos(fcy)) + 8(cos(fc2;) + COS 



i{ky)) - 4(cos(2fc:r + ky) 



- 2ky)) 



4cos(fc.,;) + 4cos(fcy) + 2(cos(2fcj: + ky 



- cos 



i{2kx ky) 



- cos{kx — 2ky) + cos{kx + 2ky)) 



+40S3(i) 



2 / ,/\ 2 



tV ft' 



t 



8(cOs(fca;) + COS 

COS(fcy)) 

12(cos(fcK) + cos(fcy)) + (2cos(2fc: 



S(fcy)) + (2 C0S(2fc^ - fcy) + 2 C0S(2fc^ + fcy) + 2 COS(fc; 



- 2 cos(fc2; - 2fcy)) + 4(cos(fc3 
■ tV ft'' 



- ky) + 



2COS(2fca; + fcy) + 2COS(fc2; + 2 fcy ) 



- 2 cos(fc2: — 2fcy)) 



(C2) 
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APPENDIX D: SPIN WAVE 
RENORMALIZATION FACTOR 

We consider the terms in the transformed spin only 
Hamiltonian of order l/S"^ as a perturbation to the lowest 
order spin wave Hamiltonian. Within this framework the 
perturbation gives quartic terms in boson creation and 
annihilation operators, for which we have to diagonalize 
the set of 2 X 2 matrices of elements (0|7k-ff47k |0), where 
(7k, 7^) e {"k, /3k}- 

We first review the situation when t' and t" are set 
equal to zero: 



SJlit) + b]b,) + 2{a,b, + alb]), 

<i,3> 

2 ^ 

<hj> 



(^alaittibj + Qibljbjbj + aja|ai 



Mblb^b, 



2a\aib]bj. 



(Dl) 

Here, we work to leading order in t/U , for which .Ji{t) = 
At^ /U . Note that, within this appendix, we use a com- 



r(4) 



pact notation: Hl'^J^^^^ H2 and Hl^^^^^ rii. ^ 
out that the 2x2 perturbation matrix is proportional to 
the identity matrix. We find: 



(4) 



Ha. It turns 



2N 
-2{u 



vl) J2 ("k'Wk' (2 cos(fc;) + 2 cos(fcp) 

k' 

UkVk(2cos{kx) + 2cos(fcy))^ 2v^, 

k' 



+4Mki'k ^ Uk'i'k' (2 cos(fc^ — fc^) + 2 cos(fcyA:y)) 

k 



(D2) 

from which the magnon energy correction Ak — Aek/ck 
can be computed. Ak converges numerically to give 

Ak ~ 0.1579, (D3) 

uniformly over the Brillouin zone, as found previously in 
Ref. [Hi. 

When t' or t" are brought into the picture the leading 
effects are to create second and third nearest neighbor 
interactions J2{t') = 4(t')V?7 and J^it") = A{t"f/U 
producing new contributions to the quartic Hamiltonian: 



+a\i 



E 

<<jj>> 



a]aiaia\ 



E 



<<<ij>>> 



a\aiaia] 



+a\a}-a]aj + a\a\aiaj + a\a]ajaj 



(D4) 



The contribution coming from t' reads: 



5el 



2N 



Mt') 2{ul + vl)Y,'^vl, 
l + vl)Y^2{ul,Zl,)x 

k' 

(2cos{{k^ + ky) - {k'^ + k'y)) + 



2cOs((fc2: 



iK -k'y)) 



+ (Wk + ^'k)E("k'+«k')x 
k' 

2(cos(fc; + fc;)+cos(fc^-/c;) 
+2(cos(A:j. + ky) + cos(/sj; — ky))x 

Y.{ul,+vl)2vl 



(D5) 



while that from t" reads: 



2{ul + vl)J2K' 

k' 

("k + vl) ^ 2{ul, + vl,)(2 cos(2fc, - 2k'J 



+ 2cos{2ky ~2ky] 

+ («k + «k) E("k' + «k')2( cos(2fc^) + cos(2fc;; 



+2(^cos(2k,,) + cos(2fey)^ (ul, + vl,)2vl, 

(D6) 

We now find that the total contribution of the terms com- 
ing from t, t' and t" is no longer independent of wave 
vector, but the two-fold degeneracy for each k remains. 
The evolution of Zc(k) over the Brillouin zone is shown 
in Fig. [2] for the parameter set (|38p . 

Dividing the contributions to Zc(k) into three parts 
coming from t, t' and t" and denoting an average over 
the Brillouin zone by (•)k we find 




0.187, 
0.019, 
0.013, 



(D7) 



using the parameter set ([38|l . Note that here AJ^. and 
A'^ do not factorize analytically in Eq. (jPSP or (jD6p . 
rather, for each wave vector, we calculate AJ^ = Se[^/e[^ 
and similarly for A'^. Note also that as t' and t" become 
non-zero, even the first term evolves. This is because 
is itself a function of t' and t" and so it changes as these 
parameters are switched on. 

As discussed in Section IHI A41 incorporating the k- 
dependence of Zc(k) makes the fitting procedure some- 
what computationally cumbersome and slow from a CPU 
speed point of view. There are two reasons for this. 
Firstly, it is due to the need to constantly redo the "inter- 
nal" sum over k' in Eqs. (|D2p . (|D5P and (|D6P whenever 
the t, t' , t" and U parameters are readjusted. Secondly, 
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150 200 250 300 350 400 



FIG. 7: (Color online) Evolution of (Afc)j, as a function of the 
number of points in the Brillouin zone 



200 250 

L 



FIG. 8; (Color online) Evolution of L x (A^)^ as a function 
of the number of points in the Brillouin zone. 



APPENDIX E: PARTIALLY CONSTRAINED 
SPIN- WAVE FITS 

For a number of compounds, such as Sr2Cu02Cl2, 
Bi2Sr2CaCu208+5 (BSCCO, Bi-2212) and 
YBa2Cu306+5 (YBCO), values of t' - -0.3 eV, 
t" ~ 0.2 eV, and ?7 > 10 eV are commonly used2^. 
However, these values are not truly universal among the 
cuprates. For example, in Ref. [2^, local density approx- 
imation (LDA) calculations do predict variations from 
one compound to the other. In particular, for La2Cu04, 
t' /t = —0.17 is obtained, and which corresponds to the 
value used in ARPES data analysis^!. 

A recent ARPES experiment by Yoshida et alr^ on 
doped La2-xSr3:Cu04 finds that t' varies slightly between 
compositions. For x = 0.03, they obtain t' /t « —0.2 and 
t" /t' « 0.5. To the best of our knowledge, there exist no 
ARPES measurements on undopcd La2Cu04, but one 
may extrapolate to t' /t = —0.21 for this composition. 
In this context, we have performed a constrained fit to 
the spin- wave dispersion data of Ref. 0], imposing t' /t = 
-0.2, t"/t' = 0.5 and t = 0.25 eV and only allowing U 
to vary. We were unable to find a reasonable value of U 
with such constraints. Indeed, wc found U ^ 1.6 cv with 
a poor quality of fit. 

We then allowed both t and U to vary. The motivation 
being that the effective ratios t' /t = —0.2 and t" /t' = 0.5 
could be well determined by ARPES, while the value 
of t determined from the electronic bandwidth could be 
strongly renormalized. The best fit then obtained was 
t ^ 0.3430 cV and U = 2.55 eV, hence U/t = 7.42. 
Compared to the fit with free t, t' , t" and U values, the 
overall x^, with 



to compound this problem the CPU time is further in- 
creased by noting that the sums over k' in Eqs. (|D2[) . 
(|D5[) and (|D6[) converge somewhat slowly with system 
size, k' in equations (jD2|D5|D6P This is illustrated in 
Fig. [7l for a magnetic Brillouin zone of linear dimension 
L, for the parameter set ([38)) (here the lattice parameter 
is taken to be unity) Re-plotting the data in Fig. [HI we 
show that (Afc) scales with system size as 



(Afe),(L) = - + (A,),(oo) , 

where a is a constant, from which we find: 
(Afe),.(oo) ~ 0.219 . 



(D8) 



(D9) 



The width of the dispersion of A(k) over the Brillouin 
zone, CT(A(k)), scales in a similar way and we eventually 
find 



^(oo) ^ 0.33x10 



-3 



(DIO) 



Finite size effects lead to negligible corrections compared 
to the least squares fitting procedure, for L ^ 100. A 
magnetic Brillouin zone of this size was used to make the 
fits described in the main text. 



(experimental uncertainty)^^ 



(El) 



where Ei is n'the magnon energy data point, is increased 
by roughly 25% for the fit with t' /t and t" /t' constrained 
compared with the fit with the values in Eq. ([55[) . 

This value of t ~ 0.34 eV agrees roughly with those 
found by Coldea et ali^, and the value U « 2.55 is in- 
between that of Ref. and the one found in this work 
and reported in Eq. ([55[) . In other words, for the con- 
strained values of t' and t" used in such a fit, the values 
are sufficiently "small" that the results for t' = t" = 
of Ref. are more or less recovered. In other words, 
the ratio t/U is a crucial parameter in the fit. Working 
with a fit with t'/t and t" /t fully constrained leads to 
a (reduced) ratio t/U compared to the unconstrained fit 
of Eq. ([38[) . Yet, this reduced U ~ 2.55 eV also gives a 
reduced t = 0.34 ev, with the result that while the un- 
constrained fit has t/U = 0.126, the constrained fit has 
t/U = 0.34/2.55 = 0.133 (I just calculated) and Coldea 
et at finds t/U ~ 0.135. Given that t'/t" is essentially 
the same in the constrained and unconstrained fits, one 
can think of the constrained fit as an intermediate point 
in switching on t' and t" at constant t'/t". Indeed, we 
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find the following progression of values, moving from the 
Coldea et al. results, the constrained fit results and the 
unconstrained fit results: U ~ 2.3 eV, 2.55 eV and 3.34 
eV, t = 0.3 eV, 0.34 eV and 0.42 eV, and t/U = 0.135, 
0.133, and 0.126. This progression of parameters ap- 
pears monotonous. Given that we expect t' and t" to 
frustrate Neel ordering, we also expect the metal insu- 



lator transition also to appear for smaller values oi t/U 
with increasing t' and t" . Our results are reassuring in 
this respect because we find the fitted t/U value getting 
progressively small as t' and t" are turned on, consistent 
with the system remaining in the insulating Neel ordered 
state. 
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Proof 



[ri,A4] = _J_ ^ \n^cTc\^Cjahj^, ol.<T'nk,a'\ (-1) 
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X/ljsO-^,a'(-l) 
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Since ( — 1)' = —(—1)-' for nearest neighbors, we get the 
final result and can generalize it for T-i and Tq. 

™ Proof : 

Same calculation than in Eq. (|48l49l50p . but for second or 
third nearest neighbors, ( — 1)' = (^l)"*- 



